Skip to content

Add mass matrix support for FIRK and MIRK - #521

Open
SebastianM-C wants to merge 2 commits into
SciML:masterfrom
SebastianM-C:smc/mass_matrix
Open

Add mass matrix support for FIRK and MIRK#521
SebastianM-C wants to merge 2 commits into
SciML:masterfrom
SebastianM-C:smc/mass_matrix

Conversation

@SebastianM-C

@SebastianM-C SebastianM-C commented Jul 2, 2026 •
edited
Loading

Copy link
Copy Markdown
Member

via unprojected collocation for index-1 DAEs based on Ascher & Spiteri (1994)

Checklist

  • Appropriate tests were added
  • Any code changes were done in a way that does not break public API
  • All documentation related to code changes were updated
  • The new code follows the
    contributor guidelines, in particular the SciML Style Guide and
    COLPRAC.
  • Any new documentation only uses public API

Additional context

Add any other context about the problem here.

@github-actions

github-actions Bot commented Jul 2, 2026 •
edited
Loading

Copy link
Copy Markdown
Contributor

Benchmark Results (Julia v1.11)

Time benchmarks
master e5ea434... master / e5ea434...
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK2() 0.666 ± 0.0067 s 0.652 ± 0.012 s 1.02 ± 0.021
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK3() 13.2 ± 0.43 ms 13.2 ± 0.27 ms 1 ± 0.039
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4() 2.82 ± 0.1 ms 2.83 ± 0.098 ms 0.997 ± 0.05
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK5() 3.42 ± 0.12 ms 3.38 ± 0.11 ms 1.01 ± 0.048
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK6() 1.64 ± 0.086 ms 1.63 ± 0.082 ms 1.01 ± 0.073
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false) 1.57 ± 0.084 ms 1.6 ± 0.097 ms 0.984 ± 0.079
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true) 3.3 ± 0.64 ms 3.22 ± 1.2 ms 1.03 ± 0.43
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false) 0.0428 ± 0.0048 s 0.0435 ± 0.0044 s 0.985 ± 0.15
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true) 0.0669 ± 0.027 s 0.0653 ± 0.0012 s 1.02 ± 0.42
Simple Pendulum/IIP/Shooting(Tsit5()) 0.187 ± 0.085 ms 0.182 ± 0.083 ms 1.03 ± 0.66
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK2() 0.824 ± 0.015 s 0.838 ± 0.013 s 0.983 ± 0.024
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK3() 18 ± 8 ms 17.8 ± 7 ms 1.01 ± 0.6
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK4() 3.53 ± 0.21 ms 3.52 ± 0.26 ms 1 ± 0.095
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5() 4.31 ± 0.31 ms 4.3 ± 0.31 ms 1 ± 0.1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK6() 2 ± 0.18 ms 2.01 ± 0.15 ms 0.998 ± 0.12
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false) 3.63 ± 1.3 ms 3.61 ± 1.2 ms 1 ± 0.5
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true) 7.23 ± 6.4 ms 6.95 ± 7 ms 1.04 ± 1.4
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false) 0.106 ± 0.0038 s 0.101 ± 0.005 s 1.06 ± 0.064
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true) 0.167 ± 0.0079 s 0.158 ± 0.0075 s 1.06 ± 0.071
Simple Pendulum/OOP/Shooting(Tsit5()) 0.655 ± 0.11 ms 0.644 ± 0.069 ms 1.02 ± 0.21
time_to_load 7.83 ± 0.13 s 7.4 ± 0.059 s 1.06 ± 0.02
Memory benchmarks
master e5ea434... master / e5ea434...
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK2() 0.389 M allocs: 0.0443 GB 0.389 M allocs: 0.0443 GB 1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK3() 0.0436 M allocs: 4.88 MB 0.0436 M allocs: 4.88 MB 1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK4() 15.7 k allocs: 1.65 MB 15.7 k allocs: 1.65 MB 1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK5() 22.1 k allocs: 2.04 MB 22.1 k allocs: 2.04 MB 1
Simple Pendulum/IIP/BoundaryValueDiffEqMIRK.MIRK6() 12.8 k allocs: 1.05 MB 12.8 k allocs: 1.05 MB 1
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = false) 25.5 k allocs: 1.83 MB 25.5 k allocs: 1.83 MB 1
Simple Pendulum/IIP/MultipleShooting(10, Tsit5; grid_coarsening = true) 0.049 M allocs: 3.39 MB 0.049 M allocs: 3.39 MB 1
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = false) 0.553 M allocs: 0.0535 GB 0.553 M allocs: 0.0535 GB 1
Simple Pendulum/IIP/MultipleShooting(100, Tsit5; grid_coarsening = true) 0.833 M allocs: 0.0778 GB 0.833 M allocs: 0.0778 GB 1
Simple Pendulum/IIP/Shooting(Tsit5()) 4.64 k allocs: 0.224 MB 4.64 k allocs: 0.224 MB 1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK2() 0.89 M allocs: 0.984 GB 0.89 M allocs: 0.984 GB 1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK3() 0.0932 M allocs: 24.8 MB 0.0932 M allocs: 24.8 MB 1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK4() 0.0324 M allocs: 3.96 MB 0.0324 M allocs: 3.96 MB 1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK5() 0.045 M allocs: 4.99 MB 0.045 M allocs: 4.99 MB 1
Simple Pendulum/OOP/BoundaryValueDiffEqMIRK.MIRK6() 25.3 k allocs: 2.17 MB 25.3 k allocs: 2.17 MB 1
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = false) 0.142 M allocs: 10.2 MB 0.142 M allocs: 10.2 MB 1
Simple Pendulum/OOP/MultipleShooting(10, Tsit5; grid_coarsening = true) 0.266 M allocs: 18.7 MB 0.266 M allocs: 18.7 MB 1
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = false) 2.52 M allocs: 0.279 GB 2.52 M allocs: 0.279 GB 1
Simple Pendulum/OOP/MultipleShooting(100, Tsit5; grid_coarsening = true) 3.82 M allocs: 0.404 GB 3.82 M allocs: 0.404 GB 1
Simple Pendulum/OOP/Shooting(Tsit5()) 0.0373 M allocs: 1.69 MB 0.0373 M allocs: 1.69 MB 1
time_to_load 0.159 k allocs: 11.2 kB 0.159 k allocs: 11.2 kB 1

@SebastianM-C
SebastianM-C marked this pull request as ready for review July 3, 2026 14:11
@ChrisRackauckas-Claude ChrisRackauckas-Claude mentioned this pull request Jul 8, 2026
Closed
SebastianM-C and others added 2 commits July 30, 2026 18:04
@SebastianM-C @claude
Add mass matrix support for FIRK and MIRK
2427637 
via unprojected collocation for index-1 DAEs based on Ascher & Spiteri (1994)

Co-authored-by: Claude Opus 4.6 <noreply@anthropic.com>
@SebastianM-C @claude
Error on adaptive solves for DAEs and add index-1 DAE tests
e5ea434 
Defect-based mesh adaptivity cannot converge for DAEs since the
collocation interpolant is inaccurate for algebraic variables, so
requesting it now throws an ArgumentError pointing users to
`adaptive = false` or BoundaryValueDiffEqAscher.jl.

Also fix the nested FIRK out-of-place collocation kernels calling the
in-place algebraic constraint helper on the user's out-of-place f.

Tests cover MIRK2-6, RadauIIa and LobattoIIIc (expanded and nested) on
a simple index-1 DAE and on example problem 1 from Ascher & Spiteri
(1994), plus the adaptivity error. LobattoIIIa/IIIb are excluded: their
tableau structure cannot determine the algebraic stage components.

Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Reviewers

No reviews

Assignees

No one assigned

Labels

None yet

Projects

None yet

Milestone

No milestone

Development

Successfully merging this pull request may close these issues.

1 participant

@SebastianM-C